Plane- Wave  Detection: 
A  Nonlinearly  Ill- Posed 
Inverse  Problem 

William  W.  Symes 
September,  1989 


TR89-12 


Report  Documentation  Page 

Form  Approved 

OMB  No.  0704-0188 

Public  reporting  burden  for  the  collection  of  information  is  estimated  to  average  1  hour  per  response,  including  the  time  for  reviewing  instructions,  searching  existing  data  sources,  gathering  and 
maintaining  the  data  needed,  and  completing  and  reviewing  the  collection  of  information.  Send  comments  regarding  this  burden  estimate  or  any  other  aspect  of  this  collection  of  information, 
including  suggestions  for  reducing  this  burden,  to  Washington  Headquarters  Services,  Directorate  for  Information  Operations  and  Reports,  1215  Jefferson  Davis  Highway,  Suite  1204,  Arlington 

VA  22202-4302.  Respondents  should  be  aware  that  notwithstanding  any  other  provision  of  law,  no  person  shall  be  subject  to  a  penalty  for  failing  to  comply  with  a  collection  of  information  if  it 
does  not  display  a  currently  valid  OMB  control  number. 

1.  REPORT  DATE 

SEP  1989 

2.  REPORT  TYPE 

3.  DATES  COVERED 

00-00-1989  to  00-00-1989 

4.  TITLE  AND  SUBTITLE 

5a.  CONTRACT  NUMBER 

Plane-Wave  Detection:  a  Nonlinearly  Ill-Posed  Inverse  Problem 

5b.  GRANT  NUMBER 

5c.  PROGRAM  ELEMENT  NUMBER 

6.  AUTHOR(S) 

5d.  PROIECT  NUMBER 

5e.  TASK  NUMBER 

5f.  WORK  UNIT  NUMBER 

7.  PERFORMING  ORGANIZATION  NAME(S)  AND  ADDRESS(ES) 

Computational  and  Applied  Mathematics  Department  ,Rice 

University, 6100  Main  Street  MS  134, Houston, TX, 77005-1892 

8.  PERFORMING  ORGANIZATION 

REPORT  NUMBER 

9.  SPONSORING/MONITORING  AGENCY  NAME(S)  AND  ADDRESS(ES) 

10.  SPONSOR/MONITOR'S  ACRONYM(S) 

11.  SPONSOR/MONITOR'S  REPORT 
NUMBER(S) 

12.  DISTRIBUTION/AVAILABILITY  STATEMENT 

Approved  for  public  release;  distribution  unlimited 

13.  SUPPLEMENTARY  NOTES 

14.  ABSTRACT 

15.  SUBIECT  TERMS 

16.  SECURITY  CLASSIFICATION  OF: 

17.  LIMITATION  OF 
ABSTRACT 

18.  NUMBER 

OF  PAGES 

32 

19a.  NAME  OF 
RESPONSIBLE  PERSON 

a.  REPORT 

unclassified 

b.  ABSTRACT 

unclassified 

c.  THIS  PAGE 

unclassified 

Standard  Form  298  (Rev.  8-98) 

Prescribed  by  ANSI  Std  Z39-18 


Plane- Wave  Detection:  a  Nonlinearly 
Ill-Posed  Inverse  Problem 


William  W.  Symes 

Department  of  Mathematical  Sciences 
Rice  University 
Houston,  Texas  77005  U.S.A. 

September  1989 


Acknowledgements.  It  is  a  pleasure  indeed  to  thank  Professor  Gary  Roach 
for  his  magnificent  hospitality  during  the  Workshop  on  Inverse  Problems  and 
Imaging,  Ross  Priory,  16-21  July  1989;  for  the  opportunity  to  present  this 
material  there;  and  for  the  invitation  to  publish  it  in  these  Proceedings.  Pro¬ 
fessor  Guy  Chavent  of  Universite  de  Paris  IX  (Dauphine)  and  the  Program 
CEREMADE  are  also  due  thanks  for  their  support  of  the  research  reported 
here,  which  was  also  partially  supported  by  Grant  No.  N00014-J-89-1 119 
and  the  National  Science  Foundation  under  grant  DMS  86-03614. 


1  Introduction 


It  is  commonplace  that  inverse  problems  of  applied  mathematics  are  often  ill- 
posed.  Such  a  problem  generally  takes  the  (idealized  )  form  of  a  functional 
equation  involving  a  map  (f>,  representing  the  physical  connection  between 
model  and  data 

<f){x)  =  z 

with  x  (the  model,  or  solution)  and  2  (the  data)  ranging  over  suitable  function 
spaces,  on  which  the  map  (j>  is  defined.  Such  a  problem  is  said  to  be  ill-posed 
when — roughly  speaking — the  solution  x  does  not  depend  continuously  on 
the  data  2.  A  common,  and  much-studied,  source  of  ill-conditioning  is  the 
appearance  of  arbitrary  small  singular  values  of  the  linearized  map  D(j)(x). 
We  might  call  the  instability  of  the  solution  resulting  from  small  singular 
values  of  D<f>  linear  ill-posedness.  See  e.g.  Tikonov  and  Arsenin  (1974). 

Some  inverse  problems  exhibit  a  form  of  effective  ill-posedness  of  quite 
different  nature.  This  second  pathology  may  appear  when  the  functional 
equation  is  re-written  as  a  best-fit  problem: 

min  ||^(x)  —  z\\ 

with  a  suitable  choice  of  norm  j|  ||.  Such  reformulation  might  be  motivated 
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by  the  recognition  that  the  range  of  (j>  does  not  contain  all  interesting  data 
points  z — -i.e.  that  the  data  is  likely  to  be  inconsistent,  perhaps  because  of 
simplifications  in  the  physics  used  to  formulate  the  map  cf>.  If  <f>  is  “sufficiently 
noninear,”  the  cost  function 

J(x-,z)  =  ||<?i(2:)  -  z\\ 

may  be  nonconvex,  with  many  local  minima  over  a  prescribed  set  of  ad¬ 
missible  models  x.  For  various  reasons — e.g.  computational  cost  of  global 
search — one  may  be  forced  to  regard  any  local  minimum  of  J  as  a  “solution.” 
Then  the  existence  of  many  (local)  minima  must  be  regarded  as  a  form  of 
instability  of  the  solution — i.e.  nonlinear  ill-posedness. 

One  inverse  problem  exhibiting  a  severe  case  of  nonlinear  ill-posedness 
is  the  velocity  estimation  problem  of  reflection  seismology.  For  an  extensive 
discussion  of  this  rather  specialized  problem,  with  many  references,  we  direct 
the  reader  to  the  monograph  by  Santosa  and  the  author  (Santosa  and  Symes, 
1989).  In  the  following  pages  we  will  mostly  discuss  instead  a  much  simpler 
problem,  the  plane  wave  detection  problem ,  which  shares  the  essential  math¬ 
ematical  features  of  the  velocity  estimation  problem  without  carrying  the 
conceptual  baggage  of  reflection  seismology.  In  the  final  section  we  will  de- 
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scribe  briefly  a  version  of  velocity  estimation,  to  make  plausible  this  sharing 
of  features. 

A  deep  understanding  of  nonlinear  ill-posedness  and  related  matters  is  to 
be  had  through  G.  Chavent’s  theory  of  quasiconvex  sets  in  Hilbert  space 
(Chavent,  1980).  The  plane-wave  detection  problem  is  treated  from  the 
point  of  view  of  Chavent’s  theory  in  Symes  (1989).  In  the  present  paper 
we  describe  just  the  basic  properties  of  a  simple  best-fit  formulation  of  the 
detection  problem  (§2),  and  then  indicate  how  the  cost  function  can  be  con- 
vexified  and  smoothed  (§3).  The  fourth  section  begins  by  describing  the 
acoustic  model  of  reflection  seismology,  then  presents  several  simplifications 
and  approximations  leading  to  a  problem  recognizably  similar  to  plane  wave 
detection.  We  give  a  brief  discussion  of  this  problem,  referring  the  interested 
reader  to  other  papers  for  more  information. 
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2  The  Output  Least-Squares  Formulation 


We  suppose  that  the  function  z(f,t)  is  a  sampling  on  {rj  =  0}  of  the  three- 
dimensional  field  U(£,r),t ): 

z(Z,t)  =  u(t,o,t),  -1<£<1,  teJR. 

We  also  suppose  that  U  is  (approximately)  a  plane- wave  moving  at  speed  1, 
i.e. 

U (£,  rj,  t )  =  u(t  —  £  sin  9  —  rj  cos  6)  +  n(£,  rj,  t) 

for  a  waveform  u(t),  an  angle  9,  and  a  noise  component  n(£,  rj,  t ).  The  plane- 
wave  detection  problem  is: 

Given  z(£,  t),  estimate  the  waveform  u(t )  and  the  angle  0  (or  equiva¬ 
lently  the  sine  s  =  sin  6). 

Because  one  believes  that  the  noise  n(£,0,t)  is  small  in  the  mean-square 
sense,  or  for  statistical  reasons  (Tarantola,  1987),  one  may  assume  that  z  E 
T2(  [— 1, 1]  x  1R),  and  ask  that  the  data  2  be  predicted  optimally  in  the  least- 
squares  sense  by  a  pair  u(t),s.  That  is,  making  also  the  a  priori  assumption 
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that  supp  u  C  [0, 1],  i.e.  u  G  L2[ 0, 1],  we  define 

#b  «](£,*)  =  u(t  -  si) 

<t>  :  [— 1)  1]  x  T2[0, 1]  — >  L2{  [—1, 1]  x  1R) 

Set 

E0  =  [-1, 1]  x  Z2[0, 1] 

F  =  L2([- 1,1]  xIR) 

with  obvious  norms.  Then  an  output-least  squares  formulation  of  the  plane- 
wave  detection  problem  is: 

min  ||#r]  -  zfF  . 

xE&o 

Remark  One  could  imagine  that  U  represents  the  far-field  signal  of  an  un¬ 
derwater  acoustic  source.  Then  the  problem  becomes  that  of  estimating  the 
waveform  emitted  by  the  source  and  its  direction  relative  to  a  line  of  receivers. 
With  this  interpretation,  our  problem  becomes  a  caricature  of  an  important 
problem  in  ocean  acoustics.  We  note  that  our  point  of  view  is  quite  different 
from  that  taken  in  the  ocean  acoustics  literature. 

An  obstacle  to  the  study  of  the  above  least-squares  problem  is  immedi¬ 
ately  evident:  As  defined,  the  map  (j>  is 

•  continuous,  but 
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•  nowhere  locally  uniformly  continuous, 


•  hence  a  fortiori  nowhere  differentiable. 
In  fact 


<j>[s  +  6s,  u](£,  t)  -  <j>[s,  u](f ,  t)  =  u(t  -  (s  +  8s)£)  -  u(t  -  s£) 

and  the  stated  properties  of  <f>  follow  from  familiar  properties  of  the  shift  on 
£2(IR).  Thus  we  cannot  study  the  dependence  of  the  solution  on  the  data 
via  the  implicit  function  theorem  or  related  tools,  nor  can  we  use  Newton’s 
method  or  its  relatives  to  compute  minima  with  any  confidence  of  conver¬ 
gence. 

The  map  <j>  becomes  of  class  C 2  if  its  domain  is  restricted,  say,  to 

E2  :=  [— 1, 1]  x  Hq[0,  1]  . 

This  restriction  does  not  cure  the  problem  of  its  delinquent  features,  of  course. 
To  begin  with,  in  any  F-neighborhood  of  any  consistent  data  point  a  =  <f>(x), 
there  exist  consistent  data  z'  =  <f>(x'),  with  as  large  as  one  likes. 

Therefore  some  regularization  of  the  optimization  problem  is  necessary,  in 
order  that  the  solution  depend  stably  on  the  data,  even  though  <j>  is  not  a 
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smoothing  operator!  This  is  a  rather  trivial  sort  of  instability,  however;  the 
actual  state  of  affairs  is  much  worse.  We  will  establish: 


There  exist  consistent  data  z  =  (f)(x o)  with  ||xo||B2  <  1  for  which  the 
problem 

min 

has  at  least  two  (local)  solutions. 

Thus  restricting  the  H2- size  of  the  solution  does  not  restore  well-posedness 
to  the  best-fit  version  of  the  detection  problem,  even  for  noise-free  data! 

Set  z(£,t)  =  u0(t)  =  ax(t)sinut  with  u  and  a  to  be  determined,  and 
X  £  Co°(0, 1)  fixed.  Then 


Hix)  -  z\\F 


z  =  <j>[x0],  x0  —  [0,  u0]  . 


For  u  €  L2[0, 1]  ,  s  €  [—1, 1]  ,  x  =  [5, it], 


(*>#*])  =  a  M  dt  u{t  —  s£)x(t)  smut 

-  •/>/  dt  u(t)x{t  +  s£) (sin  ut  cos  us£  +  cos  cot  sino;.s£) 
sinws  f 

=  2 a -  /  dt  u(t)x(t)  smut 

LOS  J 

+  J  d£  J  dt  u(t)(x{t  +  s£)  —  x(0)[s^na;^  coswsx  +  w  cos  ut  sina;s^] 
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There  is  a  uniform  estimate  for  derivatives  of  x  of  order  <  2: 


x(k\t  +  sO-x(k)(t)\<c\s 


for  5  €  [—1, 1],  £  €  [—1, 1],  t  6  IR  and  k  <  2.  Accordingly  an  integration- 
by-parts  argument  shows  that  the  second  term  is  bounded  in  absolute  value 

by 

au>  (7|s|  ||u||tf2[0(1]  • 


Thus 

{z,  <f>[x])F  =  2—J^-(u,  u0)l2[o,i]  +  0(u~2a\s\  IN||^2[0j1])  • 

So 


II^N  -  z\\f  =  +  \\z\\f  -  z)f 


>  2 


II  1 1 2  II  1 1 2  Sill  LOS  .  . 

llUo||L2[0ii]  +  ||«||l2[0,1]  _  2  us  (U’U0/L2[0,1) 


—  Cu>  2a|s 


>  s 


1  - 


ullr/2[o,i] 
sinews  I 


uis 

—  Cu>~2a\s\  II w 


u 


^[0,1] 


+  “0 


L2  [0 


.1]) 


\H2  [0,1] 


Integration-by-parts  shows  that  there  exist  C\ ,  C2  >  0  so  that 


°2  ^  +  of)  “  liU°lli-2[0,l]  -  ' 
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The  hypotheses  of  the  statement  to  be  proved  allow  us  to  require  that 
IMI/P[o,i]  <  1-  Thus  for  s0  =  %  and  x  =  [60,u]  we  get 


\\<f>[x]  -  z\\F  >  2a2  -Ccnroj-3 

>  a2Ci(2  —  e) 


provided  that 


Cl£>^+C* 


u  auJ 


Now  ||^o||#2[0)1]  <  1,  provided  |u;|  >  1  and  au 2  <  K  for  suitable  K  >  0. 
Thus  we  take  a  =  Ku~2  so  that  the  above  condition  becomes 


C\e  >  (2 C2  +  CxK-1)^-1 


and  is  satisfied  for  any  choice  of  e  >  0  as  soon  as  u  is  large  enough. 
Now  consider  the  special  choice  u  —  au0.  Then 


|| 4>[x\  —  z\\  <  2 


1  + a2  -2a 


smus 


us 


IIwoIIl2[o,i]  +  Cau  2I5Iq  IIuo||^2[0i1] 


Choose  Si  =  j1,  x  =  [si,o:uo];  then  the  above  is 


<  2 


1  +  Q2  - 


4a 


57T  J 
4a 


IKH^o,!]  +  SCaanu  3  ||m0||H2[o,i] 

C2 


<  2a2  1  +  a2  —  —  (C\  H - +  3aCairu  3  . 

07T  J  \  U  / 
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Choose  a  =  |j:;  we  get 


for  u  large  enough. 

Now  choose  e  =  (57t)-2.  Then  we  have  shown  that,  for  m  sufficiently 
large 

(i)  for  any  u  with  |M|„2[01]  <1,  s0  =  z  =  [s0,u]: 

Ml*]  ~  z\\ f  >  a2Ci(2  -  e) 

(ii)  for  S!  =  g,  U\  =  £u0,  Xi  =  [sx,u] 

life]  -  z\\f  -  q2Ci(2  -  2e)  . 

Since  any  continuous  path  from  Xi  to  x0  must  pass  over  the  set  {x  =  [s,  u] 
we  have  shown  that  the  set 

{x  =  [s,«]  :  ||fe  -  z\\F  <  a2C\(2  -  3/2e)} 

is  not  connected ;  in  particular  the  component  containing  X\  is  disjoint  from 
the  component  containing  Xq.  The  connected  component  of  Xq  is  contained 
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in 

Co  =  ja;  =  [a,  u]  :  |s|  <  —  —  <!>,  ||w||^2[0il]  <  1  j 

for  a  suitable  choice  of  8  >  0,  and  the  connected  component  of  x\  is  contained 
in 

Ci  =  =  [s,  n]  :  —  +  ^  <  s  <  1, 11^11^2^,1]  —  ^  j1 

which  follows  from  the  uniform  continuity  of  <f)  on  [—1,1]  x  if1  [0,1]  and 
the  compactness  of  the  injection  H 2  — ►  H1.  The  sets  Co  and  C\  are  closed, 
bounded,  and  convex  in  E 2,  hence  weakly  closed,  whence  follows  the  existence 
of  a  local  minimizer  in  each.  In  particular,  we  have  established  the  existence 
of  a  local  minimum  distinct  from  Xo,  as  required.  q.e.d. 

It  is  easy  to  extend  this  reasoning  to  generate  examples  with  any  number 
of  local  minima  whatsoever.  Thus  even  the  restriction  of  (j)  to  a  ball  in  E2 
does  not  suffice  to  render  the  output  least-squares  problem  well-posed  in  the 
nonlinear  sense. 
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3  The  Coherence  Reformulation 

To  motivate  the  next  step,  regard  the  data  z(£,  t)  of  the  detection  problem 
as  a  number  of  independent  time  series  measurements,  parameterized  by  the 
location  £  along  the  x-axis.  This  point  of  view  reflects  accurately  the  way 
such  measurements  are  actually  made — in  reality,  only  a  “few”  receivers  may 
be  deployed,  and  the  sample  rate  in  £  is  far  lower  than  that  in  t. 

The  plane-wave  hypothesis  implies  that  these  time  series  are  not  inde¬ 
pendent,  but  are  tied  rigidly  together  by  the  time  delay  rate  s.  The  difficulty 
described  in  the  last  section  is  also  a  consequence  of  this  rigidity.  It  is  simply 
very  difficult  to  match  all  of  the  time  series  at  once  with  any  but  precisely 
the  “right”  time  delay  rate.  Any  other  choice  results  in  large  mismatch 
somewhere. 

Our  solution  to  the  quandary  of  the  hyperactive  behaviour  of  the  model  is 
to  relax  it.  We  allow  independent  models  of  the  various  traces,  constrained  by 
a  penalty  for  deviation  from  the  plane-wave  hypothesis.  The  penalty  is  not 
“capital  punishment,”  as  in  the  least-squares  approach;  models  inconsistent 
with  the  plane-wave  hypothesis  are  permitted,  but  required  to  pay  a  “fine” 
related  to  their  deviation.  The  “fine  schedule,”  ie.  penalty  weight,  is  a  very 
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important  determinant  of  efficiency  of  such  a  scheme. 


Precisely,  we  introduce  a  cover  <j>  of  the  map  <j>,  that  is  a  diagram  of  maps 
and  spaces 

4> 


E 


a 


<t> 


P 


E 


which  partially  commutes: 

if  4>(x)  €  'R.(P),  then 
p(x)  =  p(<f>(a(x)))  . 

We  also  assume  that  ft  is  injective.  As  we  shall  see,  it  is  possible  for  such  a 
partially  commuting  diagram  to  be  constructed  with  a  differentiable  <j&,  even 
though  <f>  is  not  differentiable. 

First  note  the  consequence:  if 

5  =  Hz) 

and  <f>(x)  =  z,  then  x  =  a(x)  satisfies  </>(x)  =  z,  since  (3  is  injective.  That 
is,  the  functional  equation  for  <f>  has  the  “same”  solutions  as  the  functional 
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equation  for  <f>,  provided  that  the  data  for  the  former  is  the  image  under  /? 
of  the  data  for  the  latter. 

The  space  E ,  as  a  set,  is 


E  =  [-l,l]x{fi€/f’([-l,l]x[-l,3])  : 

=  0,  i  <  -1C| 

du 

=  o,  <>i  +  iei} . 

We  shall  actually  need  a  family  of  Hilbert  spaces  Ea,  a  >  0,  consisting  of  the 
set  E  equipped  with  the  norm 


5,  u 


2 

E,cr 


du 


dt 


+  <?  \sr  + 


du 


dt 


(here  the  subscript  “0”  denotes,  as  before,  the  norm  in  L2(  [-1, 1]  x  [-1, 3] ).) 
We  shall  also  need  a  family  of  maps 


k  :  K  ->  F  :=  F  ®  F 


defined  by 

k[s,u] 

Finally,  a  and  (3  are  given  by 


a[s,  u\ 
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m  = 


(A 


Obviously  /?  is  injective.  Also,  if  4>[s,u\  6  71(13),  then 


du  du 
d£~Sdt 


=  0 


so  | q(£,t)  =  u(t  +  s()  for  some  u  £  £2(IR). 

Because  of  the  assumptions  about  the  support  of  u,  supp  u  C  [0,1]. 
Moreover 


a[s,  u]  =  [3,  u]  . 


Necessarily 

]>[x,u]  =  /3(<j)[s,u})  . 


As  noted  above,  for  consistent  data  z  £  7Z(f3),  the  functional  equation  for 
<j)  has  the  “same”  solution  as  that  for  <f>.  On  the  other  hand,  <j>c  :  Ea  — >  F 
is  differentiable — even  real- analytic.  Of  course  the  diagram  above  could  not 
commute  if  all  maps  but  <j>  were  differentiable.  In  fact,  a  is  not  C1.  However, 
application  of  a  may  be  viewed  as  a  “postprocess,”  to  get  from  <f>~l(z )  to 
4>~1(z),  and  does  not  entail  solving  a  functional  equation. 

Because  <j)a  is  C 2  or  better,  local  well-posedness  of  the  least-squares  prob- 
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lem 

min  U(7[s,u]  -  ill  —Ja[s,u-,z] 

(s,u)  I  II 

for  near-consistent  z  follows  from  a  coercivity  estimate  on  the  derivative  D(j). 
Thus  assume  that 

<M-So,u0]  =  S0  =  /3(z0)  . 

We  shall  prove  that 

\\i^,8u}\\^  <  L.[uo]  ||Z>^[s0,m][^,^u]||^ 

for  some  La  >  0.  It  follows  immediately  from  this  estimate  and  the  implicit 
function  theorem  that  Ja  has  a  unique  minimum  near  [s0,fio]  for  z  near  i0- 
The  main  step  in  the  proof  of  the  coercivity  estimate  is  the  estimate 

+  K1\\DW[s0,u0}[8s,8u]\\l 

with  Ko,Ki  independent  of  [so,fio]-  Here  W  is  the  second  component  of  <j>\\ 

and  it  is  assumed  that  ^[s0,  Wo]  €  7£(/?),  i-e.  in  particular  that  W[so?ffo]  =  0. 
Set 


du0 

2 

\Ss\2  <  Ko 

d6u 

dt 

0 

dt 

D  :  =  DW[so,  «o][<5s,  <5u] 

duo  c  d6u  d6u 

=  -mSs--dTSo  +  W 
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Pick  £i,£2  €  [—1, 1],  ^2  =  ^1  +  -5o(^2  —  £i).  Then  integrate  both  sides  of  the 
previous  equation  along  the  line  segment 

£  £  [£i5&]  —  6)) 


to  get 


du0 

dt 


(6^i)^(6  —  £i)  —  f  d£D(£i,ti  +  s0(£  —  £i))  + 

•Ni 

+  8u(^\,ti)  —  ^(^2>  ti  +  ^0(^2  —  £1)) 
According  to  the  definition  of  E,  6u  =  0  for  t  <  |£|,  so  ( t  <  3) 


|£u(£,t)|2  <  4  J  dt' 


d6u 


dt 


(£/) 


so 


(6-6)2l^12 


du0 


dt 


(6,<0 


<  3  f2di  \D{tth  +  sit  -  h))?  + 


r3  , 

dSu . 

2 

d8u  2 

Ldt 

dt 

+ 

a  «*.")  ] 

Now  integrate  both  sides  in  £21  ^1  and  £1  in  that  order  to  get 


>'2 


du0 


dt 


<12||A>||2  +  48 


ddu 


dt 


whence  we  can  take  Kp  =  72,  AT  =  18.  Thus  the  “main  step”  estimate  is 
established. 
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To  finish  the  proof  of  the  coercivity  estimate,  note  that 


dSu 

2 

<  3 

dt 

0 

P>llo  +  N2 


du0 

2 

■  1  1 9 

dSu 

2' 

dt 

+  |so| 

0 

dt 

0 

whence 
2 


du0 


dt 


[<5s,<5n]  \\\a  = 


duQ 


dt 


dSu 


dt 


+  cr2  |^|2  + 


< 


du0 


dt 


dSu 


dt 


+  72a 2  +  219a 


du0 


dt 


d8u 


dt 


+  a2  18  +  57 


du0 


dt 


\\D\\l 


Thus  the  coercivity  estimate  is  proved  with 


L<t[soi  ft0]  =  max  <  1  +  219cr2  +  72a2 


du0 


dt 


-2 


,  a2  57  +  18 


du0 


dt 


— 2  > 


As  mentioned  above,  this  estimate  ensures  the  local  existence,  uniqueness, 
and  continuous  dependence  of  the  global  minimizer  of  Ja  for  near-consistent 
data,  i.e.  data  near  the  intersection  of  the  range  of  </v  and  the  range  of  (3. 
This  local  well-posedness  result  is  just  the  beginning  of  the  story,  of  course. 
Two  questions  of  immediate  “practical”  interest  are: 


(i)  How  does  one  ensure  that  an  initial  estimate  of  x  =  [s,u]  lies  suf¬ 
ficiently  close  to  the  global  minimum  of  Ja  to  permit  convergence 
of  a  Newton-like  iterative  scheme? 
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(ii)  How  can  one  be  sure  that  the  data  is,  or  is  not,  close  enough  to 
consistent  data  that  the  local  well-posedness  result  holds? 

The  first  question  is  addressed  in  the  paper  (Symes,  1989)  using  Chavent’s 
quasiconvexity  theory.  There  we  show  that  a  suitable  initial  estimate,  and 
so  a  convergent  quasi-Newton  sequence,  may  be  constructed,  provided  that 

(a)  the  noise  level  is  sufficiently  low; 

(b)  <r  is  chosen  sufficiently  small  initially,  and  later  increased  to  provide 
the  maximum  level  of  stability. 

All  of  these  results  rely  on  estimates  which  are  doubtless  overly  conservative. 
A  major  open  problem  is  the  derivation  of  algorithms  to  estimate  appropriate 
values  of  a  and  of  the  noise  level,  hence  answering  question  (ii). 
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4  The  Seismic  Reflection  Inverse  Problem 


In  this  section  we  give  a  rather  terse  description  of  the  inverse  problem  of 
reflection  seismology,  and  explain  its  relation  to  the  plane  wave  detection 
problem  discussed  in  the  preceding  sections. 

Reflection  seismology  is  an  active  remote  sensing  technique.  It  currently 
yields  the  most  highly  detailed  images  available  for  study  of  the  earth’s  crust 
to  depths  of  15-30  km.  The  method  was  developed  by  the  petroleum  in¬ 
dustry,  but  enjoys  some  use  in  academic  geophysics  (tectonics,  continental 
margins)  and  engineering  geophysics  (structure  placement,  groundwater)  as 
well. 

The  data  of  reflection  seismology  are  measurements  of  the  transient  re¬ 
sponse  of  the  crust  to  controlled  energy  sources — explosive  charges,  com¬ 
pressed  air  devices,  dropped  weights,  hydraulic  vibrators,  etc.  The  mea¬ 
surements  of  ground  velocity  (on  land)  or  pressure  (at  sea)  are  made  by  a 
(typically  linear)  array  of  sensors.  A  typical  marine  seismic  cable  used  in 
oil  and  gas  exploration  may  be  3  km.  long  and  contain  several  hundred 
hydrophones  (immersible  microphones),  wired  together  in  groups  to  produce 
signals  on  48,  96,  or  even  more  channels.  The  sample  rate  for  recording  these 
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signals  is  250-1000  samples  per  second,  and  the  typical  digitally  recorded  sig¬ 
nal  may  be  six  seconds  long.  The  source  may  be  repeatedly  stimulated,  each 
10  seconds  or  so,  as  the  exploration  vessel  steams  along  towing  the  cable.  A 
brief  calculation  shows  that  the  amount  of  data  recorded  in  a  day’s  work  by 
such  a  ship  is  truly  phenomenal.  In  fact,  the  geophysical  exploration  industry 
has  been  for  many  years  the  single  largest  industrial  user  of  digital  magnetic 
tape.  An  excellent  and  easily  accessible  overview  of  reflection  seismology  and 
other  geophysical  prospecting  methods  is  provided  by  Dobrin  (1975). 

The  extraction  of  useful  information  from  this  vast  quantity  of  data  is  par¬ 
ticularly  challenging  because  of  the  complex  physics  of  seismic  wave  propaga¬ 
tion.  In  fact,  this  physics  is  only  poorly  understood.  It  is  generally  accepted 
that  some  modification  of  linear  elasticity  might  be  suitable,  as  the  limited 
time-scale  and  bandwidth  of  seismic  waves  conspire  to  hide  or  average  the 
microscopic  mechanics  of  rocks.  In  fact  most  of  the  research  literature,  and 
almost  all  production  software,  is  based  on  the  acoustic  model,  which  repre¬ 
sents  the  earth  by  a  fluid.  The  local  linear  response  of  a  fluid  is  characterized 
by  a  density  p(x )  and  sound  velocity  c(x),  both  functions  of  spatial  position. 
If  the  energy  source  is  assumed  isotropic  (it  isn’t),  and  small  (on  the  scale  of 
the  dominant  wavelength,  say)  and  positioned  at  xs  then  the  excess  pressure 
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p(x,t )  obeys 


1 


d2p(X,t)  -  v  •  -L-vp(x,t)  =  nm  x  -  x.) 


p(x)c(x)2  dt 2  v  ’  '  p(x) 

p(x,  t)  =  0  ,  t  <<  0  . 


In  principle,  a  boundary  condition  is  also  needed  at  the  earth/air  interface 
(x3  =  0:  the  earth’s  surface  is  flat  on  the  scale  of  the  exploration  seismic 
experiment!),  but  for  simplicity  we  shall  assume  that  p,  c  are  constant,  at 
fixed,  known  values  for  X3  <  0  (“up”  in  geophysics!). 

This  initial  value  problem  represents  the  simplest  model  of  seimic  wave 
propagation  retaining  most  of  its  notable  features.  In  terms  of  the  model,  the 
basic  problem  of  reflection  seismology  is  to  recover  the  mechanical  descriptors 
of  the  earth,  i.e.  p(x)  and  c(x),  from  the  idealized  surface  seismogram 


s[p,c](xi,x2,t)  =  p(xi,x2,0,t)  . 


To  exhibit  the  close  relation  with  the  plane-wave  detection  problem,  three 
further  simplifications  are  useful.  Each  limits  the  utility  of  the  model  further 
—  though  the  end  result  can  still  be  used  successfully  to  process  field  seismic 
data! 

First,  we  assume  “the  earth  is  layered,”  i.e.  that  p  and  c  depend  only  on 
z  =  x3.  This  assumption  is  a  reasonable  approximation  as  sediments  which 
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later  become  rock  are  laid  down  in  horizontal  layers,  and  retain  a  large  degree 
of  lateral  homogeneity  in  many  areas.  Because  of  the  assumed  translation 
invariance,  we  can  replace  the  seismogram  by  a  transform  in  the  horizontal 
variables.  We  use  the  Radon  transform,  suitably  scaled: 


P(f,z,t)  —  J  J  dxxdx2  p(xi,x2,z,t  +  £xi) 


which  obeys 

z’ ()  =  mHz  ~z‘)  =  0 

P(£,z,t)  =  0,  t  «  0  . 

The  plane-wave  seismogram  is 


=  PU,0,t)  . 


Note  that  the  I.V.P.  for  P  is  hyperbolic  where  |£|c  (*)  <  L 

The  second  simplification  is  based  on  the  recognition  that  the  density  and 
velocity  in  the  sedimentary  earth  exhibit  a  dichotomy  of  scales :  small-scale 
fluctuation  in  the  meters-tens  of  meters  range,  corresponding  to  seasonal  and 
short-term  changes  in  weather  patterns  and  depositional  environment,  and 
long  scale  fluctuations  in  the  hundreds  or  thousands  of  meters  corresponding 
to  geological,  epochal  changes.  It  was  recognized  early  that  the  short-scale 
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fluctuations  could  be  viewed  usefully  as  perturbations  about  the  long-scale 
trends  (though  the  mathematical  justification  for  this  linearization  is  recent 
and  still  partial).  Accordingly  write 

p  —  ps  +  pr 

c  —  cs  -j-  cr 

where  the  subscripts  “s”  and  “r”  stand  for  “smooth”  and  “rough”  respec¬ 
tively.  Regular  perturbation  about  cs,  ps  gives 

DS\ps",  cs][/>r,  cr]  =  SP |2=0 

where 

1  /  _  2\  d2SP  _d_}_  dSP_  1  2ct  d2P  d  pr  dP 
ps  \cs  )  dt 2  dz  ps  dz  ps  dt 2  dz  pi  dz 

and 

SP  =  0,  <«0. 

Note  that  SP  is  linear  in  pT  and  cr,  but  nonlinear  in  ps  and  cs. 

The  third  simplification  is  to  use  geometric  optics  to  approximate  the 
solution  of  the  above  perturbational  boundary  value  problem.  This  step  is 
justified  by  the  recognition  that  the  wavelengths  measured  in  seismic  records 
are  quite  short  compared  to  scale  of  the  slowly-varying  coefficients  (cs  and 
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ps ).  After  some  work  (e.g.  Santosa  and  Symes,  1988)  one  obtains 


DS\pSj  c5] [pr ,  cr] (£?  <£)  ~  A(^ps^  cSj  £) 


'>r+-+^2'’' 


/M 


°C(6*) 


Here  £  is  the  inverse  two-way  traveltime ,  defined  implicitly  by 


rC(€,<)  /  1 

t  =  2  I  dzj- — -£2 


=  2/ 


\{z) 


and  A(ps,  cs;  £,  i)  is  an  overall  amplitude  factor,  determined  by  ray  geometry. 

Close  examination  of  the  last  two  equations  reveals  the  following  compar¬ 
isons  with  the  plane-wave  detection  problem: 


(1)  The  Radon  parameter  (“slowness”)  £  plays  the  role  of  receiver  co¬ 
ordinate  £;  t  is  time  in  both  cases. 

(2)  The  smooth  background  velocity  cs  plays  the  role  of  the  direction 
sine  s,  in  that  it  determines  a  £- dependent  change  of  variables:  for 
the  reflection  seismology  problem, 


for  the  plane-wave  detection  problem,  the  time  shift  t  h- ►  t  +  s£. 

(3)  pr / ps  and  cTjcs  play  the  role  of  the  plane  waveform  u(t). 
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In  a  number  of  recent  papers,  we  have  made  deductions  for  the  reflection 
seismology  model  similar  to  those  sketched  in  sections  1-3  for  the  plane- 
wave  detection  problem.  We  have  shown  that  the  map  (ps,cs,  pr,cr)  i— » 
DS[ps,  cs][pr,  cr]  is  either  non-differentiable  or  has  non-coercive  derivative, 
depending  on  choice  of  metric  in  the  domain.  We  have  also  posed  a  coherency 
method,  by  replacing 


Mfl  |  *(*) 

Ps{z )  CS(Z ) 

Pr(z) 

Ps(z) 


rd(t,0 


and  requiring  that  the  reflectivities  ri  and  satisfy  the  coherency  conditions: 

|noC  =  ° 
d 

% r-o<=o 


These  are  analogous  to  the  condition 


U(U)  =  0 


of  the  plane-wave  detection  problem.  The  fit-to-data  and  coherency  con¬ 
straints  can  be  combined  into  a  least-squares  problem,  just  as  was  shown  in 
Section  3.  We  have  built  a  Fortran  code  to  solve  the  resulting  least-squares 
problem. 
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The  theory  of  this  approach  to  the  reflection  seismology  problem,  includ¬ 
ing  well-posedness  estimates,  extension  to  the  original  acoustic  model  (rather 
than  the  simplified,  linearized  plane  wave  model),  and  preliminary  numer¬ 
ical  experiments,  are  reported  in  Symes  (1988a, b).  Application  of  the  co¬ 
herency  technique  to  field  seismic  data  is  accomplished  in  Symes  and  Caraz- 
zone  (1989). 
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